###reduced version
use.formula = as.formula(paste('in.out.diff.dictator ~ ',master.formula,sep = ''))
use.formula.reduced = as.formula(paste('in.out.diff.dictator ~ ',master.formula.reduced,sep = ''))

use.vars = c('in.out.diff.dictator',master.use.vars)

#UO
use.dat.uo = dat.uo[,use.vars]
use.dat.uo = use.dat.uo[complete.cases(use.dat.uo),]

##N per city for weighted regression below
city.N = tapply(use.dat.uo$income,use.dat.uo$index.city,length)
city.N = as.data.frame(city.N)
city.N$index.city = row.names(city.N)
use.dat.uo = merge(use.dat.uo,
                    city.N)

use.dat.uo$index.city.numeric = as.numeric(use.dat.uo$index.city)

reg.uo.dict.reduced = lm(use.formula.reduced, data = use.dat.uo)
clustered.se = cl(use.dat.uo, reg.uo.dict.reduced, use.dat.uo$index.city)
reg.uo.dict.reduced$se = clustered.se
reg.uo.dict.reduced.weighted = lm(use.formula.reduced, data = use.dat.uo,
                                  weights = city.N)
brl.se <- brl(use.formula.reduced, dataframe=use.dat.uo, 
              clusterName="index.city.numeric", method="brl")[,2]
reg.uo.dict.reduced.brl = reg.uo.dict.reduced
reg.uo.dict.reduced.brl$se = brl.se


reg.uo.dict = lm(use.formula, data = use.dat.uo)
clustered.se = cl(use.dat.uo, reg.uo.dict, use.dat.uo$index.city)
reg.uo.dict$se = clustered.se  ##note that this does not show up in the object until apsrtable() is applied so summary(reg.uo.dict) will show hte wrong standard errors
reg.uo.dict.weighted = lm(use.formula, data = use.dat.uo,
                           weights = city.N)
brl.se <- brl(use.formula, dataframe=use.dat.uo, 
              clusterName="index.city.numeric", method="brl")[,2]
reg.uo.dict.brl = reg.uo.dict
reg.uo.dict.brl$se = brl.se


#Non UO
use.dat.sec = dat.non.uo[,use.vars]
use.dat.sec = use.dat.sec[complete.cases(use.dat.sec),]

##N per city for weighted regression below
city.N = tapply(use.dat.sec$income,use.dat.sec$index.city,length)
city.N = as.data.frame(city.N)
city.N$index.city = row.names(city.N)
use.dat.sec = merge(use.dat.sec,
                    city.N)

use.dat.sec$index.city.numeric = as.numeric(use.dat.sec$index.city)

reg.sec.dict.reduced = lm(use.formula.reduced, data = use.dat.sec)
clustered.se = cl(use.dat.sec, reg.sec.dict.reduced, use.dat.sec$index.city)
reg.sec.dict.reduced$se = clustered.se
reg.sec.dict.reduced.weighted = lm(use.formula.reduced, data = use.dat.sec,
                                   weights = city.N)
brl.se <- brl(use.formula.reduced, dataframe=use.dat.sec, 
                               clusterName="index.city.numeric", method="brl")[,2]
reg.sec.dict.reduced.brl = reg.sec.dict.reduced
reg.sec.dict.reduced.brl$se = brl.se

reg.sec.dict = lm(use.formula, data = use.dat.sec)
clustered.se = cl(use.dat.sec, reg.sec.dict, use.dat.sec$index.city)
reg.sec.dict$se = clustered.se
reg.sec.dict.weighted = lm(use.formula, data = use.dat.sec,
                           weights = city.N)
brl.se <- brl(use.formula, dataframe=use.dat.sec, 
              clusterName="index.city.numeric", method="brl")[,2]
reg.sec.dict.brl = reg.sec.dict
reg.sec.dict.brl$se = brl.se


###############
outtable = apsrtable(reg.uo.dict.reduced,
                     reg.uo.dict,
                     reg.sec.dict.reduced,
                     reg.sec.dict, 
                     Sweave = T,
                     coef.names = coef.names,
                     model.names = c('UO','UO','Non UO','Non UO'),
                     notes = ''
                    )
writeLines(
  outtable, 'output/appendix/Table_A7.tex')



outtable = apsrtable(reg.uo.dict.reduced.weighted,
                     reg.uo.dict.weighted,
                     reg.sec.dict.reduced.weighted,
                     reg.sec.dict.weighted, 
                     Sweave = T,
                     coef.names = coef.names,
                     model.names = c('UO','UO','Non UO','Non UO'),
                     notes = ''
                    )
writeLines(
  outtable, 'output/appendix/Table_A16.tex')

outtable = apsrtable(reg.uo.dict.reduced.brl,
                     reg.uo.dict.brl,
                     reg.sec.dict.reduced.brl,
                     reg.sec.dict.brl, 
                     Sweave = T,
                     coef.names = coef.names,
                     model.names = c('UO','UO','Non UO','Non UO'),
                     notes = ''
              )
writeLines(
  outtable, 'output/appendix/Table_A19.tex')



###zelig for effects of segregation

##UO
reg.uo = zelig(use.formula, data = use.dat.uo,
		model = 'ls',
		cite = F)

	
#expected values
outputs.uo = matrix(nrow = length(dissim.scale.uo), ncol = 3)
for(i in dissim.scale.uo){
	row.i = which(dissim.scale.uo == i)
	outx = setx(reg.uo, outgroup.dissim.yeshiva = i,  outgroup.yeshiva = uo.homo.value,
		jerusalem = 0,  	
    demo1.sex = sex, 
		demo1.age = age, 
		demo1.ethnicity = ethnicity,
		demo2.left_right = ideology,
		income = income.value, 
		college = college.value, 
		nonjewish_pcnt = nonjewish_pcnt.value,
		immigrant = 0
	)
	
	sims = sim(reg.uo,outx, sims = n.sims)
	outs = sims$qi[[1]]
	outputs.uo[row.i,1] = mean(outs)
	outputs.uo[row.i,2] = quantile(outs,low.conf)
	outputs.uo[row.i,3] = quantile(outs,high.conf)
	}

	
###Non UO
reg.sec = zelig(use.formula, data = use.dat.sec,
		model = 'ls',
		cite = F)		
		
outputs.sec = matrix(nrow = length(dissim.scale.sec), ncol = 3)
for(i in dissim.scale.sec){
	row.i = which(dissim.scale.sec == i)
	outx = setx(reg.sec, outgroup.dissim.yeshiva = i, outgroup.yeshiva = sec.homo.value,
		jerusalem = 0,
		demo1.sex = sex, 
		demo1.age = age, 
		demo1.ethnicity = ethnicity,
		demo2.left_right = ideology,
		income = income.value, 
		college = college.value, 
		nonjewish_pcnt = nonjewish_pcnt.value,
		immigrant = 0
	)
	
	sims = sim(reg.sec,outx, sims = n.sims)
	out = summary(sims)$stats["Expected Values: E(Y|X)"]
	outs = sims$qi[[1]]
	outputs.sec[row.i,1] = mean(outs)
	outputs.sec[row.i,2] = quantile(outs,low.conf)
	outputs.sec[row.i,3] = quantile(outs,high.conf)
	}


	
###UO	
pdf('output/Figure_2_regarding_seg_uo.pdf',
    width = pdf.width,
    height = pdf.height)    
	
par(mar = c(4,4.5,.25,.25))
par(las = 1)
plot(0,0,
	ylim = dict.ylims,
	xlim = xlims.dissim.uo,
	type = 'n',
	xlab = 'Segregation',
	ylab = 'Ingroup Bias',
	cex.axis = cex.axis,
	cex.lab = cex.lab)
	
abline(h = 0,
		col = 'gray',
		lty = 2)
lines(dissim.scale.uo,
	outputs.uo[,1],
	lwd = 4,
	col = line.color
	)
lines(dissim.scale.uo,
	outputs.uo[,2],
	lty = 3,
	lwd = 3,
	col = line.color
	)
lines(dissim.scale.uo,
	outputs.uo[,3],
	lty = 3,
	lwd = 3,
	col = line.color
	)
	
dev.off()
	
##Non UO
pdf('output/Figure_2_regarding_seg_non_uo.pdf',
    width = pdf.width,
    height = pdf.height)    

#par(mar = c(4.5,5.5,.25,.25))
par(mar = c(3,4.5,.25,.25))

par(las = 1)	
plot(0,0,
	ylim = dict.ylims,
	xlim = xlims.dissim.sec,
	type = 'n',
  xlab = '',
  ylab = 'Ingroup Bias',
	cex.axis = cex.axis,
	cex.lab = cex.lab)
	
abline(h = 0,
		col = 'gray',
		lty = 2)
lines(dissim.scale.sec,
	outputs.sec[,1],
	lwd = 4,
	col = line.color
	)
lines(dissim.scale.sec,
	outputs.sec[,2],
	lty = 3,
	lwd = 3,
	col = line.color
	)
lines(dissim.scale.sec,
	outputs.sec[,3],
	lty = 3,
	lwd = 3,
	col = line.color
	)
	
dev.off()


#####################
###zelig for effects of homogeneity
##UO
outputs.uo = matrix(nrow = length(homo.scale.uo), ncol = 3)
for(i in homo.scale.uo){
	row.i = which(homo.scale.uo == i)
	outx = setx(reg.uo,  outgroup.yeshiva = i, outgroup.dissim.yeshiva = mean.seg.value,
		jerusalem = 0,
		demo1.sex = sex, 
		demo1.age = age, 
		demo1.ethnicity = ethnicity,
		demo2.left_right = ideology,
		income = income.value, 
		college = college.value, 
		nonjewish_pcnt = nonjewish_pcnt.value,
		immigrant = 0
	)
	
	sims = sim(reg.uo,outx, sims = n.sims)
	out = summary(sims)$stats["Expected Values: E(Y|X)"]
	outs = sims$qi[[1]]
	outputs.uo[row.i,1] = mean(outs)
	outputs.uo[row.i,2] = quantile(outs,low.conf)
	outputs.uo[row.i,3] = quantile(outs,high.conf)
	}


outputs.sec = matrix(nrow = length(homo.scale.sec), ncol = 3)
for(i in homo.scale.sec){
	row.i = which(homo.scale.sec == i)
	outx = setx(reg.uo, outgroup.yeshiva = i, outgroup.dissim.yeshiva = mean.seg.value,
		jerusalem = 0,
		demo1.sex = sex, 
		demo1.age = age, 
		demo1.ethnicity = ethnicity,
		demo2.left_right = ideology,
		income = income.value, 
		college = college.value, 
		nonjewish_pcnt = nonjewish_pcnt.value,
    immigrant = 0
	)
	
	sims = sim(reg.sec,outx, sims = n.sims)
	out = summary(sims)$stats["Expected Values: E(Y|X)"]
	outs = sims$qi[[1]]
	outputs.sec[row.i,1] = mean(outs)
	outputs.sec[row.i,2] = quantile(outs,low.conf)
	outputs.sec[row.i,3] = quantile(outs,high.conf)
	}

	
#UO	
pdf('output/Figure_2_regarding_proportion_uo.pdf',
    width = pdf.width,
    height = pdf.height)    
	
par(mar = c(4.5,4.5,.25,.25))
par(las = 1)	
plot(0,0,
	ylim = dict.ylims,
	xlim = xlims.homo.uo,
	type = 'n',
	xlab = 'Outgroup Proportion',
	ylab = 'Ingroup Bias',
	cex.axis = cex.axis,
	cex.lab = cex.lab)

abline(h = 0,
		col = 'gray',
		lty = 2)
lines(homo.scale.uo,
	outputs.uo[,1],
	lwd = 4,
	col = line.color
	)
lines(homo.scale.uo,
	outputs.uo[,2],
	lty = 3,
	lwd = 3,
	col = line.color
	)
lines(homo.scale.uo,
	outputs.uo[,3],
	lty = 3,
	lwd = 3,
	col = line.color
	)

dev.off()
	
##Non UO
pdf('output/Figure_2_regarding_proportion_non_uo.pdf',
    width = pdf.width,
    height = pdf.height)    
	
par(mar = c(3,4.5,.25,.25))

par(las = 1)
plot(0,0,
	ylim = dict.ylims,
	xlim = xlims.homo.sec,
	type = 'n',
  xlab = '',
	ylab = 'Ingroup Bias',
	cex.axis = cex.axis,
	cex.lab = cex.lab)

abline(h = 0,
		col = 'gray',
		lty = 2)


lines(homo.scale.sec,
	outputs.sec[,1],
	lwd = 4,
	col = line.color
	)
lines(homo.scale.sec,
	outputs.sec[,2],
	lty = 3,
	lwd = 3,
	col = line.color
	)
lines(homo.scale.sec,
	outputs.sec[,3],
	lty = 3,
	lwd = 3,
	col = line.color
	)
	
dev.off()

